Hands on Exercise 2b- Second Order Spatial Point Patterns Analysis Methods

2b 2nd Order Spatial Point Patterns Analysis Methods

Built By Deepika Raja

1. Getting Started

1.a. Learning Outcome

Second-order spatial point pattern analysis examines the spatial relationships between points in a pattern, specifically focusing on how the presence of one point influences the location of others. It goes beyond simply describing the overall density of points (first-order effects) by investigating clustering, dispersion, or randomness at various spatial scales.

2. Data

• Child Care Services data from data.gov.sg, a point feature data providing both location and attribute information of childcare centres.

• Master Plan 2019 Subzone Boundary (No Sea), a polygon feature data providing information of URA 2019 Master Plan Planning Subzone boundary data.

3. Installing and Loading the R packages

pacman::p_load(sf, terra, spatstat, 
               tmap, rvest, tidyverse)

4. Data Import and Prepration

mpsz_sf <- st_read("data/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>% st_transform(crs = 3414)

Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source 
  `C:\deepikarcodes\ISSS626-GAA\Hands_on_Exercises\Hands_on_Ex02\data\MasterPlan2019SubzoneBoundaryNoSeaKML.kml' 
  using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

extract_kml_field <- function(html_text, field_name) {
  if (is.na(html_text) || html_text == "") return(NA_character_)
  page <- read_html(html_text)
  rows <- page %>% html_elements("tr")
  value <- rows %>%
    keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
    html_element("td") %>%
    html_text2()
  if (length(value) == 0) NA_character_ else value
}

mpsz_sf <- mpsz_sf %>%
  mutate(
    REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
    PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
    SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
    SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
  ) %>%
  select(-Name, -Description) %>%
  relocate(geometry, .after = last_col())

mpsz_cl <- mpsz_sf %>%
  filter(SUBZONE_N != "SOUTHERN GROUP",
         PLN_AREA_N != "WESTERN ISLANDS",
         PLN_AREA_N != "NORTH-EASTERN ISLANDS")

write_rds(mpsz_cl, "data/mpsz_cl.rds")

childcare_sf <- st_read("data/ChildCareServices.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = 3414)

Reading layer `CHILDCARE' from data source 
  `C:\deepikarcodes\ISSS626-GAA\Hands_on_Exercises\Hands_on_Ex02\data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

st_crs(mpsz_cl) == st_crs(childcare_sf)

[1] TRUE

library(sf)
library(ggplot2)
library(readr)

# Load cleaned subzone and childcare datasets
mpsz_cl <- read_rds("data/mpsz_cl.rds")
childcare_sf <- st_read("data/ChildCareServices.kml") %>%
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = st_crs(mpsz_cl))  # Ensure same CRS

Reading layer `CHILDCARE' from data source 
  `C:\deepikarcodes\ISSS626-GAA\Hands_on_Exercises\Hands_on_Ex02\data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

# Create the map
ggplot() +
  geom_sf(data = mpsz_cl, fill = "gray", color = "black", size = 0.2) +
  geom_sf(data = childcare_sf, color = "black", size = 1.2, alpha = 0.7) +
  labs(title = "Childcare Services Across Singapore Subzones",
       subtitle = "All layers aligned to EPSG:3414 (SVY21)",
       caption = "Source: URA Master Plan & Data.gov.sg") +
  theme_minimal()

5. Interactive mapping with tmap()

tmap_mode("view")

tm_shape(childcare_sf) +
  tm_dots()

6. Switching back to plot

tmap_mode("plot")

childcare_ppp <- as.ppp(childcare_sf)

class(childcare_ppp)

[1] "ppp"

summary(childcare_ppp)

Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 11 decimal places

Mark variables: Name, Description
Summary:
     Name           Description       
 Length:1925        Length:1925       
 Class :character   Class :character  
 Mode  :character   Mode  :character  

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units

sg_owin <- as.owin(mpsz_cl)

class(sg_owin)

[1] "owin"

plot(sg_owin)

7. Let’s extract the childcare events that are located within Singapore:

childcareSG_ppp = childcare_ppp[sg_owin]

childcareSG_ppp

Marked planar point pattern: 1925 points
Mark variables: Name, Description 
window: polygonal boundary
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units

pg <- mpsz_cl %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_cl %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_cl %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_cl %>%
  filter(PLN_AREA_N == "JURONG WEST")

par(mfrow=c(2,2))
plot(st_geometry(pg), main = "Ponggol")
plot(st_geometry(tm), main = "Tampines")
plot(st_geometry(ck), main = "Choa Chu Kang")
plot(st_geometry(jw), main = "Jurong West")

pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)

childcare_pg_ppp = childcare_ppp[pg_owin]
childcare_tm_ppp = childcare_ppp[tm_owin]
childcare_ck_ppp = childcare_ppp[ck_owin]
childcare_jw_ppp = childcare_ppp[jw_owin]

childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")

par(mfrow=c(2,2))
plot(unmark(childcare_pg_ppp.km), 
  main="Punggol")
plot(unmark(childcare_tm_ppp.km), 
  main="Tampines")
plot(unmark(childcare_ck_ppp.km), 
  main="Choa Chu Kang")
plot(unmark(childcare_jw_ppp.km), 
  main="Jurong West")

8. Analysing Spatial Point Process Using G-Function

The G function measures the distribution of the distances from an arbitrary event to its nearest event. In this section, you will learn how to compute G-function estimation by using Gest() of spatstat package. You will also learn how to perform monta carlo simulation test using envelope() of spatstat package

8.a Choa Chu Kang planning area

Computing G-function estimation:

set.seed(1234)

Compute G-function using Gest() of spatat package:

G_CK = Gest(childcare_ck_ppp, correction = "border")
plot(G_CK, xlim=c(0,500))

8.b. Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1 = The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Monte Carlo test with G-fucntion:

G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(G_CK.csr)

8.c Tampines planning area

Computing G-function estimation:

G_tm = Gest(childcare_tm_ppp, correction = "best")
plot(G_tm)

8.d Complete Spatial Randomness (CSR) Test

To verify the spatial pattern of childcare services in Tampines, a hypothesis test is carried out:

• Null Hypothesis (H₀): Childcare services in Tampines are randomly distributed.

• Alternative Hypothesis (H₁): Childcare services in Tampines are not randomly distributed.

The null hypothesis will be rejected if the p-value < 0.001 (significance level).

The following code chunk executes this hypothesis test:

G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(G_tm.csr)

9. Analyzing Spatial Point Patterns with the F-Function

The F-function (or empty space function, F(r)) estimates the likelihood that a randomly chosen location is within distance r of the nearest point in a spatial pattern. It can also be interpreted through its hazard rate h(r).

This analysis uses the Fest() function from the spatstat package to estimate F(r), and envelope() for conducting a Monte Carlo simulation test to assess spatial randomness.

9.a. Choa Chu Kang Planning Area

To compute the F-function for Choa Chu Kang, the following code uses Fest() from the spatstat package.

F_CK = Fest(childcare_ck_ppp)
plot(F_CK)

9.b. Testing for Complete Spatial Randomness
To validate the observed spatial patterns, a hypothesis test is carried out. The hypotheses are defined as:

• H₀: Childcare services in Choa Chu Kang are randomly distributed.

• H₁: Childcare services in Choa Chu Kang are not randomly distributed.

  The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The test is performed using the Monte Carlo method with the F-function:

F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(F_CK.csr)

9.c. Tampines planning area

Computing F-function estimation Monte Carlo test with F-fucntion:

F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)

9.d. Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns, a hypothesis test is conducted. The hypotheses are stated as follows:

• H₀: Childcare services in Tampines are randomly distributed.

• H₁: Childcare services in Tampines are not randomly distributed.

The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The code chunk below is used to perform the hypothesis testing.

F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.

plot(F_tm.csr)

10. Analysing Spatial Point Process Using K-Function

The K-function measures the number of events found within a specified distance from any given event.

10.a Choa Chu Kang Planning Area

Computing K-function estimate

K_ck = Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")

10.b. Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns, a hypothesis test is conducted. The hypotheses are defined as follows:

• H₀: Childcare services in Choa Chu Kang are randomly distributed.

• H₁: Childcare services in Choa Chu Kang are not randomly distributed.

The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The code chunk below is used to perform the hypothesis testing:

K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")

10.c. Tampines planning area

Computing K-fucntion estimation:

K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r, 
     ylab= "K(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns, a hypothesis test is conducted. The hypotheses are defined as follows:

• H₀: Childcare services in Tampines are randomly distributed.

• H₁: Childcare services in Tampines are not randomly distributed.

The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The code chunk below is used to perform the hypothesis testing:

K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(K_tm.csr, . - r ~ r, 
     xlab="d", ylab="K(d)-r", xlim=c(0,500))

11. Analysing Spatial Point Process Using L-Function

The L-function is an alternative transformation of the K-function, making it easier to interpret spatial clustering or dispersion.

11.a. Choa Chu Kang Planning Area

Computing L-function estimation:

L_ck = Lest(childcare_ck_ppp, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

11.b. Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns, a hypothesis test is conducted. The hypotheses are stated as follows:

• H₀: Childcare services in Choa Chu Kang are randomly distributed.

• H₁: Childcare services in Choa Chu Kang are not randomly distributed.

The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The code chunk below is used to perform the hypothesis testing:

L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

11.c. Tampines planning area

Computing L-fucntion estimate:

L_tm = Lest(childcare_tm_ppp, correction = "Ripley")
plot(L_tm, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

11.d. Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns, a hypothesis test is conducted. The hypotheses are stated as follows:

• H₀: Childcare services in Tampines are randomly distributed.

• H₁: Childcare services in Tampines are not randomly distributed.

The null hypothesis (H₀) will be rejected if the p-value < α (0.001).

The code chunk below will be used to perform the hypothesis testing:

L_tm.csr <- envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(L_tm.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500))